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Abstract One goal of helioseismology is to determine the subsurface structure 
of sunspots. In order to do so, it is important to understand first the near-surface 
effects of sunspots on solar waves, which are dominant. Here we construct sim- 
plified, cylindrically-symmetric sunspot models, which are designed to capture 
the magnetic and thermodynamics effects coming from about 500 km below 
the quiet- Sun T5000 — 1 level to the lower chromosphere. We use a combi- 
nation of existing semi-empirical models of sunspot thermodynamic structure 
(density, temperature, pressure): the umbral model of Maltby et al. (1986) and 
the penumbral model of Ding and Fang (1989). The OPAL equation of state 
tables are used to derive the sound speed profile. We smoothly merge the near- 
surface properties to the quiet-Sun values about 1 Mm below the surface. The 
umbral and penumbral radii are free parameters. The magnetic field is added to 
the thermodynamic structure, without requiring magnetostatic equilibrium. The 
vertical component of the magnetic field is assumed to have a Gaussian horizontal 
profile, with a maximum surface field strength fixed by surface observations. The 
full magnetic field vector is solenoidal and determined by the on-axis vertical 
field, which, at the surface, is chosen such that the field inclination is 45° at the 
umbral-penumbral boundary. We construct a particular sunspot model based on 
SOHO/MDI observations of the sunspot in active region NOAA 9787. The helio- 
seismic signature of the model sunspot is studied using numerical simulations of 
the propagation of f, pi, and p2 wave packets. These simulations are compared 
against cross-covariances of the observed wave field. We find that the sunspot 
model gives a helioseismic signature that is similar to the observations. 

Keywords: Sun: Helioseismology, Sun: sunspots. Sun: magnetic fields 



1. Introduction 

The subsurface structure of sunspots is poorly known. Previous attempts to 
use helioseismology to determine the subsurface properties of sunspots have 
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mainly been done under the assumptions that the sunspot is non-magnetic (it is 
usuaUy treated as an equivalent sound-speed perturbation) and that it is a weak 
perturbation to the quiet-Sun. Neither of these two assumptions is justifiable. 
The helioseismolog y results have been — perhaps unsurprisingly — confusing and 
contradictory (e.g., Gizon et al. . 20091 ). The effects of near-surface magnetic and 



structural perturbations on solar waves are strong is more easily dealt with 
in numerical simulations. Examples of su ch simulations of wave propagatio n 

through prescribed sunspot models include Cameron. Gizon. and Duvalll (l2008l) : 

Hanasogg ( 2008f ): [Khomenko. Collados. and Felipe! ( ~008f ) : Parchevskv and Kosovichev^ 



(|2009f ). Necessary for all these attempts are appropriate model sunspots. There 
are numerous sunspot models available, some of which have be en used i n helio - 

seismic studies. For recen t reviews of sunspot models see, e.g., ISolanki' (2004 

Thomas and Weiss ( 2008 ), and Mprad i et al. (2009). Various authors , e.g., K homenko and Gollados^ 



(|2008l ). lMoradi and Callvl (|2008[ ). and lMoradi. Hanasoge. and Gallvl (j2009[) . have 
constructed magnetohydrostatic par ametric sunspot mo dels for use in helioseis- 
mology. In a previous paper ( Cameron. Gizon. and Duvall . 200&). we considered 
a self-similar magnetohydrostatic model. 

Although the deep structure of sunspots is of the utmost interest, it is likely 
to be swamped in the helioseismic observation s unless we are able to accuratel y 
model and remove the surface effects (see, e.g., Gizon. Birch, and Spruitl . l2010[) . 
The aim of this paper is to construct a simple parametric sunspot model, which 
captures the main effects of the near-surface layers of sunspots on the waves. 

In principle, the 3D properties (magnetic field, pressure, density, temperature, 
Wilson depression, flows) of sunspots can be i nferred in the photo sphere and 
above using spectropolarimetric inversions (e.g.. iMathew et~al. . I2OO3I) . However, 
today, these inversions are unavailable for most sunspots. In such circumstances, 
we choose to construct sunspot models that are based upon semi-empirical 
models of the vertical structure of typical sunspots. In producing our models, we 
treat separately the thermodynamic structure and the magnetic field: we do not 
require magnetostatic equilibrium. There are several reasons for not requiring 
the model to be hydrostatic. The first is that we are much more interested in the 
waves then in whether the background model is magnetodydrostatic. For this 
reason we are more interested in geting for example, the sound speed, density, 
and magnetic field to be close to those which are observed. A second reason 
is that we are not including non-axisymettric structure in our model, which 
certainly affects the force-balance. A third reason is that sunspots have both 
Evershed and moat flows, the existence of which implies a net force and so 
indicates that the sunspot is not strictly magnetohydrostatic. 

In this paper, we take the umbral and penumbral models to match those 
of existing semi-empirical models. In the absence of horizontal magnetic field 
measurements we assume that the field inclination at the umbra/penumbra 
boundary is approximately 45°. We do not consider the effects of the Evershed 
or moat fiows, although they could be included in the framework we set out. The 
surface model of the sunspot needs to be smoothly connected to the quiet-Sun 
model below the surface, and details will be provided in Section 2. Since we are 
not including the chromosphere in our simulations, we have chosen to smoothly 
transition the sunspot model above the surface to the quiet-Sun model. 
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The description for constructing models including the surface structure of 
sunspots, which will be fleshed out in Section 2, is intended to be generic. For 
illustrative purposes we choose parameter values of the sunspot model that are 
appropriate for t he sunspot in NOA A 9787, which was observed by SOHO/MDI 



in January 2002 ( Gizon et al . 20091) . In Section 3, we describe the setup of the 



numerical simulation of the propagation of f, pi, and P2 wave packets through 
the m odel sunspot. We use the SLiM code ([Cameron. Gizon. and DaifFallahl . 



20071 ). We briefly compare the simulations and the SOHO/MDI observations in 



Section 4. In Section 5, we conclude that this simple sunspot model, which is 
intended to be a good description of the sunspot's surface properties, leads to a 
good helioseismic signature and provides a testbed for future studies. 



2. Model sunspot 

The sunspot models that we describe below are cylindrically symmetric and thus 
we use a cylindrical-polar coordinate system to describe the spot, where z is the 
height and r is the horizontal distance from the axis of the sunspot. 

2.1. Thermodynamical aspects 



For this paper we use the umbral Mode l-E of Maltbv et al. ( 19861 ) and the 
penumbral model of Ding and Fane ( 19891 ). These models give the pressure and 
density as functions of height. For the quiet-Sun background in wh ich we embed 



19961 ). Zero height 



the sunspot, we use Model S (|Ghristensen-Dalsgaard et al 
{z = 0) is defined as in Model S. 

There is some freedom in choosing the depths of the umbral (Wilson) and 
penumbral depressions. We choose a Wilson depression of 550 km. More precisely, 
we place the tsooo = 1 surface of the umbra at a height 



550 km, 



(1) 



which is defined with respect to the h eight of the rsooo height of the quiet Sun 
reference model of Maltbv et al. ( 19861 ) at with zq being approximately —70 km. 
This value of 550 km for the Wilson depression produces a match between the 
density of the Maltby model and that of the quiet Sun approximately 100 km 
below the T5000 = 1 surface (see Figure [1]). 

For the penumbra, we place the T5000 = 1 surface at a height 



= zo — 150 km. 



(2) 



A penumbral d epression of 150 km i s consistent with spectropolarimetric mea- 



surements (e.g.. iMathew et aLl. 120041). 



19861) is plot ted in Figure [I ] for z 
736 km. The penumbral model of Ding and Fand (jl989l) 



The umbral model of iMaltbv et al 
z-a — 116 km = 
plotted in Figure [T] for z > Zp - 220 km = -440 km. 
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Figure 1. The red-dashed curves show the semi-empirical umbral model from'Mal tbv et al\ 
l ll986f) as a function of height. The black curves are the quiet-Sun values from model-S. The 
blue curves show the pressure the umbral model described in this paper. The panels show: 
presssure (top-left panel), density (top-right), temperature (bottom-left), and sound-speed 
(bottom-right). The vertical dashed lines are the heights where T5000 = 1 in the quiet-Sun 
(z = 20) and the umbra {z = Zu)- 



2.2. Geometrical aspects 



The umbral and penumbral models discussed above need to be smoothly embed- 
ded in Model S. In order to do so, we modify the semi-empirical models above 
and below the tsooo = 1 levels of the umbra and the penumbra, as follows. 

In our model, the pressure along the axis of the sunspot, denoted by P^, is 
such that 

InPu(z) =w;(z)lnPMaitby(2:) + (l-w(z))lnPq,(z). (3) 

In thi s expression, PMaitby(-z) is the umbral pressure given by iMaltbv et ail 
i 19861 ) for z > —736 km, and is extrapolated for z < —736 km assuming constant 
pressure scale height. The quantity Pqs{z) is the quiet-Sun pressure from Model 
S. The function w{z) is a weighting function given by 



z < Zu - 380 km, 

cos[7r <^"~^° }'^~^'> ]/2 + 1/2 Zu - 380 km < z < Zu - 80 km, 
^ 300 km J / / u _ u , 

w{z) = <( 1 Zu - 80 km < z < Zu + 20 km, (4) 

cos[7r%±2okm^i /2 _^ 1/2 z + 20 km < z < 900 km, 
^ (900 km-z„) J/ / u — ' 







900 km < z. 
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Figure 2. The red -dashed curves show the semi-empirical penumbral model from 
iDing and Fand | |1989|) . The black curves are the quiet-Sun values from Model S. The blue 
curves show the penumbral model described in this paper. Shown are, the pressure (top-left 
panel), density (top-right), temperature (bottom-left), and sound-speed (bottom-right). The 
vertical dashed lines are the heights where T5000 = 1 in the quiet-Sun (2 = zq) and the 
penumbra (2 = Zp). 



The resul ting umbral pressure is plotted in Figure [T] The pressure near z = Zu 
is that of Maltbv et al. ( 19861 ). Below, it smoothly merges with the quiet Sun 
pressure at z = Zu — 380 km = —1000 km. Above, the pressure tends to the 
quiet-Sun value and there is a significant departure from Maltbv et al. i 1986f ) 
for z > —200 km (by more than a factor of two). 

The vertical profile of the density is treated in the same manner as the pres- 
sure. The Maltby model does not explicitly contain all the properties that we 
require, so we use the OPAL tables to derive the sound speed from the pressure 
and the density. 

The penumbral pressure and the density of Ding and Fane ( 19891) are embed- 
ded in Model S in a similar way as was done for the umbra, except that the 
weighting function 'w(z) is replaced by 



w{z) 







z < Zp - 280 km. 



r (£j^_8okm_£)w2 + 1/2 Zp - 280 km < z < Zp - 80 km, 
^ 200 km J / / p — p ' 



1 



Zp - 80 km < z < Zp + 20 km, 



cos[7r%±2^^^i^]/2 -I- 1/2 Zp -I- 20 km < z < 800 km 
^ (800 km-zp) ii / P - 



800 km < z. 

The resulting penumbral pressure and density are plotted in Figure [H 



(5) 
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Figure 3. Relative weights used to combine the umbral (red), penumbral (black), and 
quiet-Sun (blue) models in the horizontal radial direction. 



Thus far we have described umbral, penumbral, and quiet-Sun models. We 
use them to form a three-dimensional cylindrically-symmetric sunspot model. 

Since we are modeling only the very near-surface layers of the sunspot (top 
1 Mm), we do not pay excessive attention to the fanning of the interfaces with 
depth. In our 3D sunspot model, each thermodynamic quantity is the product 
of a function of z and a function of r. 

The radii of the umbra and the penumbra are denoted by and rp. We 
combine the three model components (umbra, penumbra, quiet-Sun) in the r 
coordinate using the weight functions shown in Figured The transitions between 
the 1-D atmospheres have a width of 6 Mm and are described by raised cosines. 

For the sunspot in Active Region 9787, we take = 10 Mm and rp = 20 Mm 
for the umbral and penumbral radii. 

The sound speed was again reconstructed using the density, pressure, and 
the OPAL tables to ensure consistency. We note that the temperature does 
not appear in the equations and is used only for the purpose of calculating 
the properties of the spectral lines used in helioseismology (e.g., the MDI Ni 
676.8 nm line). 

2.3. Observables 



Doppler velocity is the primary observable in helioseismology. In this section, 
which is parenthetical, we consider the question of how to compute a quan- 
tity resembling the observatio ns from the suii spot model. This topic has previ- 
ously be en considered by, e.g.. IWachteii (|2008f ). Here we use the STOPRO code 
(|Solanki fl987.') ,0 synthesize, assuming LTE, the formation of the Ni 676.8 nm 
line, which is the line used by the SOHO/MDI instrument. 

We obtain a continuum intensity image near the Ni line, shown in Figure SI 
together with observations. The bright rings at the transitions between the 
umbra, penumbra, and quiet Sun are undesirable artifacts due to the simpli- 
fied transitions between the different model components. This may have to be 
addressed in a future study by adjusting the weights from Figure |3l 

W e compute the response functions for vertical velocity perturbations ( Beckers and Milk"evl .| 



19751 ) as functions of height, horizontal position, and wavelength, A. We integrate 



the response function over the wavelength bandpass [—100 mA, —34 mA] mea- 
sured from line center. This wavelength range was chosen because it is where the 



SOLA: spot_model_sp_4.tex; 2 March 2010; 9:09; p. 6 



Constructing semi-empirical sunspot models for hclioseismology 




Figure 4. Left panel: Model of the sunspot continuum intensity image near the Ni 676.8 nm 
line calculated using STOPRO. The red arcs have radii 10 Mm, 20 Mm, and 30 Mm. Noticeable 
are bright rings at the transitions between the umbra, penumbra, and quiet Sun. Right panel: 
Observed SOHO/MDI intensity image of the sunspot in AR 9787 averaged over 20 - 28 January 
2002. 
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Figure 5. The weights (being the normalized response function of the Ni 676.8 nm line) which 
allow us to create synthetic helioseismic Doppler velocity maps. The height is with reference 
to the Tsooo = 1 surface in the quiet-Sun. 



slope of the line is maximal, i.e. where the response to velocity perturbations 
is largest. The response function is normalized so that its z-integral is unity at 
each horizontal position. These are the weights (Figure [SJ to be multiplied by the 
vertical velocity and then integrated in the vertical direction to give a Doppler 
velocity map. 

Whether the bright rings affect the Doppler velocity in practice is unclear. 
We believe, however, that the response function computed here will be useful 
in studies investigating the diagnostic potential of the helioseismic signal within 
the sunspot. 

2.4. Magnetic field 

When observations of the vector magnetic field are available, they can and should 
be included in constraining the sunspot model. In this paper, we consider a 
more generic case and assume a Gaussian dependence of the vertical magnetic 
field in the radial direction. The important parameters are then the maximum 
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magnetic field at the surface, Bq, the lialf-widtli at lialf maximum (HWHM) of 
the Gaussian profile at the surface, ho, and the inclination of the magnetic field 
at the umbral-penumbral boundary. 

Explicitly, the dependence of Bz on r and z is given by 

Bzir,z) = Bzir = 0, z) exp[~{\n2y /hiz)^] (6) 

where h{z) is the half width at half maximum (HWHM) of the Gaussian. Since 
the magnetic flux at each height must be conserved, we have 

h{z) = ho^Bo/Bz{r^O,z), (7) 

where ho is the surface value of h. For the sunspot in Active Region 9787 we 
have chosen i?o = 3 kG and ho = 10 Mm. 

The function Bz{r — 0,z) is unknown. It is a goal of helioseismology to 
constrain B^ along the sunspot axis. Here we choose a two-parameter function 
as follows: 

Bz{r = 0, z) = B„e-(^/")i(-2/")i(^/"2)^ (g) 

where L{z) = 1/(1 + e^^) is the logistic function. Since dBz/dz{r = 0, z = 0) = 
— i?o/(4a), the parameter a controls the vertical gradient of B^ near the surface. 
Together with the condition V • B = 0, dBz/dz{r = 0,z = 0) determines the 
full surface vector magnetic field. We chose a = 1.25 Mm so that at z = the 
inclination of the magnetic field is 45° at the umbra/penumbra boundary, in 
agreement with observations. The parameter a2 controls the field strength at 
depth; we chose a2 = 18.4 Mm. 

Figure [B] summarizes the magnetic properties of the sunspot model. The mag- 
netic field strength increases rapidly with depth as a consequence of the choice 
of a large a2 in Equation Equation © implies in turn that the radius of the 
tube shrinks fast, from 10 Mm at the surface down to 1.3 Mm at z = — 25 Mm 
(just spatially resolved in the numerical simulations of Section 3). This sunspot 
model — as monolithic as imaginable — is only one possible choice among many. 
We note that it would be straightforward to consider model sunspots that fan 
out as a function of depth, e.g., by decreasing the parameter a2- Figure IH] also 
shows the ratio of the fast- mode speed, Vi, to the sound speed, c, along the 
sunspot axis, which is nearly unity for depths greater than 1 Mm, but increases 
very rapidly near the surface. For example, we have Uf/c = 1.9 at z = z„ and 
Vf/c — 3.9 at z = —400 km. The height at which the Alfven and sound speeds 
are equal (a = c) is z = —750 km. 

We comment that the construction of both the thermodynamic and magnetic 
properties of this model aims to get the surface properties correct. Additional 
information about the surface, when available, could easily be incorporated 
into our sunspot model. Where information is missing the assumption that the 
sunspot has similar properties to those of other sunspots is possibly the best 
that can be done. A summary of the choices made to model the sunspot in AR 
9787, which was observed by MDI, is given in Table 1. 
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Figure 6. The upper left panel shows the vertical magnetic field strength along the axis of 
the spot. The upper right panel shows the ratio of the fast mode speed to the sound speed 
also along the sunspot axis. The lower panels show the inclination (left) and strength (right) 
of the magnetic field at z = 

Table 1. Basic choices made in constructing the sunspot model for AR 9787. 



Umbral model 

Penumbral model 

Umbral radius 

Penumbral radius 

Umbral depression (rsooo = 1) 

Penumbral depression (rsooo = 1) 

On-axis surface vertical field 

Radial profile of Bz 

Field inclination at umbra/penumbra 



Maltbv et al. (1986) 

Ding and Fang (1989) 

ru = 10 Mm 

rp = 20 Mm 

20 — 2:u = 550 km 

^0 — 2p = 150 km 

Bo = 3 kG 

Gaussian with HWHM=10 Mm 
45° 



3. Numerical simulation of the propagation of waves through the 
sunspot model 



We want to simulate the propagation of planar wave packets throu gh the model 

sunsp ot an d surrounding quite-Sun using the SL iM code described in Cameron. Gizon. and Daiffallah^ 

(|2007l ) and lCameron. Gizon. and Duvalll (|2008l ). The size of the box is 145.77 Mm 

in both horizontal coordinates. We used 200 Fourier modes in each of the two 

horizontal directions; almost all of the wave energy resides in the 35 longest 

horizontal wavelengths. The vertical extent of the box is —25 Mm < z < 

2.5 Mm, where the physical domain, —20 Mm < z < 0.5 Mm, is sandwiched 

between two "sponge" layers that reduce strongly wave reflection from the 
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boundaries (jCameron. Gizon. and DuvaUl . 120081 ). We use 1098 grid points and a 
finite difference scheme in tlie vertical direction. 

Thus far the quiet-Sun reference model has been Model S which is convec- 
tively unstable. For the purpose of computing the propagation of solar waves 
through our sunspot models using linear numerical simulations, we need a stable 
background model in which to embed the sunspot model. We use Convectively 
Stabilized Model B (CSM_B) from ISchunker et al\ (|2010l ). We keep the relative 
perturbations in pressure, density, and sound speed fixed. For example, this 
means that the new pressure is given by 



P{r,z)-P{r.z)^^^f^. (9) 



where z) is the pressure as discussed above. The same treatment is applied 
to density, p, and sound speed, c. 

Separate simulations are done for f, pi, and p2 wave packets. In each case, a 
wave packet is made of approximately 30 Fourier modes, hi t — Q all the modes 
are in phase at a; = —43.73 Mm, while the sunspot is centered at x = y = 0. 
The initial conditions at t = are such that the wave packet propagates in the 
positive X direction, towards the sunspot. 

The Alfven velocity increases strongly above the photospheric layers of the 
spot (see Figure IH]). This is a problem for two reasons. The first and most 
important of these is that it causes the wavelengths of both the fast-mode and 
Alfven waves to become very large. The upper boundary, which is located at 
z = 2.5 Mm, becomes only a fraction of a wavelength away from the surface. 
This makes the simulation sensitive to the artificial upper boundary. The second 
reason is that our code is explicit and hence subject to a CFL condition: High 
wave speeds require very small time steps and correspondingly large amounts of 
computer power. In order to address this problem, we notice that we expect all 
solar waves that reach the region where the Alfven speed is high to continue to 
propagate upward out of the domain. We therefore reduce the Alfven speed in 
this region in the simulation to increase the time the waves spend in the sponge 
layer. This reduces the influence of the upper boundary and allows us to have 
a larger time step. In practice, we multiply the Lorentz force by 8c^/(8c^ -I- a?), 
where a(r, z) is the Alfven speed, and c{z) is the sound speed of the quiet Sun 
at the same height. This limits the fast-mode speed to be a maximum of three 
times the quiet-Sun sound speed at the same geometric height. Since the sound 
speed near the t — 1 layer in the spot is less than half of that of the quiet Sun 
at the same height, the maximum value of the fast mode speed in this critical 
region is approximately six times the local sound speed. The functional form 
chosen to limit the Lorentz force modifies a at and below the c — a level by less 
than 2%. This corresponds to a change in field strength of less than 30 G or a 
minimal change in the height of the c = a surface. 
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4. Preliminary comparison of simulations with observed MDI 
cross-covariances of the Doppler velocity 

There exist excellent observations (SOHO/ MDI Dopple r veloc ity) of helioseis- 
mic waves around the sunspot in AR 9787 ( Gizon et a/.l . l2009l ). The wave field 



around the sunspot can be characterized by the temporal cro ss-covariance of 



the o bserved random wave field. It has been argued (|Gizon. Birch, and Spruit 



20ld . and references therein) that the observed cross-covariance is closely re- 



lated to the Green's function (the response of the Sun to a localized source). 
Hence the observed cross-covariance is comparable to the surface vertical ve- 
lo city from initial-value numerical sim ulations of the wave packet propagation. 
In lCam eron. Gizon. and Duvalll (|2008|) we studied the propagation of an f-mode 



wave packet through a simplified magnetohydrostatic sunspot model. We found 
that we could constrain the surface magnetic field by comparing the observed 
f-mode cross-covariance with a simulated f-mode wave packet. 

Here we assess the seismic signature of the semi-empirical sunspot model 
described in Section 2 by comparing the simulation s and observations. The cross 



covari ance is constructed in the same way as in ICameron. Gizon. and Duvall 



(|2008l ). to which we refer the reader. In short, it is computed according to 



T 



C{x,y,t)^ cl>{t')cb{x,y,t' + t)dt', (10) 



where T = 7 days is the observation time, (jj is the observed Doppler velocity, 4> 
is the average of cj) over the line x = —43.73 Mm, and t is the time lag. We select 
three different wave packets (f, pi, and P2) by filtering (j) along particular mode 
ridges. As said in Section 3, for the numerical simulations we consider plane wave 
packets starting at x = —43.73 Mm and propagating in the -\-x direction towards 
the model sunspot (at the origin). The initial conditions of the simulation are 
chosen such that, in the far field, the simulated vertical velocity has the same 
temporal power as the observed cross-covariance. 

Figure [7| shows the observed cross-covariances and simulated wave packets for 
the pi modes at four consecutive time lags. The main features seen in the cross- 
covariances, the speedup of the waves across the sunspot as well as their loss of 
energy at short wavelengths, are also seen in the simulations. The details of the 
perturbed wa veforms can only be unde rstood in the context of finite- wavelength 
scattering (Gi zon. Hanasoee. and Birc h. 2006). 

For a more detailed comparison we have concentrated on the two spatial 
locations A — {x,y) — (60,0) Mm and B = (60,60) Mm, marked with crosses 
in Figure [T] The first point lies behind the sunspot, where the effects of the 
sunspot are easily noticeable. The second point is away from the scattered field 
and serves as a quiet-Sun reference. The corresponding simulated and observed 
wave packets are plotted as functions of time in Figure |S1 The match again looks 
qualitatively good. 

To proceed further we Fourier analyze the time series at points A and B: we 
consider the power spectra and phases. The temporal power spectra are shown 
in the left panels of FigurelHl The power spectrum of the simulations is similar to 
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Figure 7. (Left panels) Simulated vertical velocity of pi wave packets at times t = 0, 140, 
144, and 148 min (from top to bottom). White shades are for positive values, dark shades 
for negative values. The gray scale covers the full range of values in each row. (Right panels) 
Observed cross-covariance at the corresponding time-lags. The t = frame from the simula- 
tion is the prescribed initial condition. The red and blue crosses (hereafter points A and B 
respectively) indicate two particular locations used in the subsequent analysis (Figures |8] and 
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Figure 8. Time series from the simulations (left panels) or from the observed cross-covariance 
(right panels) for the (top to bottom) f, pi and p2 wave packets. The red curve corresponds 
to point A (the red cross) in Figure [T] the blue curve to point B (the blue cross). 

that of cross-covariance for the pi modes, less so for the f and p2 modes. We note 
that the differences between the simulated and observed power spectra are seen 
at both spatial points. We also analyzed the phase difference between the wave 
packets at the two points (A minus B) as a function of frequency (Figure H]). 
The phase shifts introduced by the sunspot are reasonably well reproduced by 
the simulations, although some differences exist especially for the f modes. 

In both the observations and the simulations we see that the power in the f 
and pi wavepackets peaks substantially below 3 mHz. The reason for this can be 
understood by noting that wave attenuation increases sharply with frequency. 
Even though the initial power spectra have most of their energy near 3 mHz, 
very little of the power at and above 3 mHz reaches the points A and B. The 
remaining low frequency waves have more sensitivity to the deeper layers and are 
thus the most important in determining the subsurface structure of the sunspot. 
It is also apparent that the sunspot more strongly 'damps' higher frequency 



SOLA: spot_model_sp_4.tex; 2 March 2010; 9:09; p. 13 



Cameron et al. 






Figure 9. Frequency analysis of the curves shown in Figure 8. The left panels shows the 
power spectra of the time series from Figure 8, with the solid lines indicating the values for the 
numerical simulations and the dots indicating the values from the observed cross-covariances. 
The red curves are for point B (away from the sunspot) abd the blue curves for point A 
(behind the sunspot). The top panels are for f mode wave packets, the middle for pi, and 
the bottom for p2. The right panels show the phase differences between points A and B, i.e. 
sunspot minus quiet-Sun. The solid lines are for the numerical simulations and the dots are 
for the observations. The error bars are root-mean-square deviations computed over 0.3-mHz 
intervals. 



waves. This is understandable since these waves have more of their energy in 
the near-surface region, where mode conversion (from the f and p modes into 
Alfven and slow-magnetoacoustic modes that propagate away from the surface) 
is efficient. Over the frequency range studied, t he pow e r 'absorption' coefficient 
incre ases with frequency (in accordance with BraunL 1995t Braun and BirchL 



20081 ). The same trend is also seen in the phase shifts, with low frequency 



waves being almost unaffected by the sunspot whilst high frequency waves have 
phase shifts which range up to 300 degrees (for the frequencies shown). A strong 
fre quency dependence of th e travel-time perturbations had alre ady been reported 
bv lBraun and Birch ( 2008 ) for ridges pi through p4 (see also lBraun and BirchL 
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l2006tlBirch et adl2009[) . The reason why the phase shifts are strongest for high- 
frequency waves is again that they have more energy near the surface. The 
possibihty of such strong phase shifts should be born in mind when interpreting 
(and measuring) hehoseismic waves near sunspots. 



5. Conclusion 

We have outhned a method for constructing magnetic models of the near-surface 
layers of sunspots by including available observational constraints and semi- 
empirical models of the umbra and penumbra. Applying this type of model to 
the sunspot of AR 9787, we showed that the observed hehoseismic signature of 
the sunspot model is reasonably well captured. Our approach was to compare 
numerical simulations of wave propagation through the model sunspot and the 
observed SOHO /MDI cross-covariances of the Doppler velocity. 

Possible improvements of the simulation include a better treatment of wave 
attenuation (esp. in plage), improved initial conditions, and inclusion of the moat 
flow. It should also be noted that improved observations of the surface vector 
magnetic field by SDO/HMI will help tune the sunspot models. 

The dominant influence of the sunspot's surface layers on hehoseismic waves 
means that it is necessary to model it very well in order to extract informa- 
tion from the much weaker signature of the sunspot's subsurface structure. The 
sunspot model used here accounts for most of the hehoseismic signal, although 
there are substantial differences that remain to be explored. In any case, the 
model provides a testbed which is sufficiently similar to a real sunspot to be 
used for numerous future studies in sunspot seismology. 

Finally, we remark that our numerical simulations of linear waves could po- 
tentially be used to interpret other hehoseismic observations than the cross- 

covariance (cf. hehoseismic holography or Fourier-Hankel an alysis), as well as 

other magnetic phenomena (e.g., small magnetic flux tubes, see Duvall. Birch, and Gizonl .[ 



20061) . 
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